#[inline(always)]
pub const fn is_prime(n: u64) -> bool {
    const NUM_BASES: usize = 11;
    #[rustfmt::skip]
    const WITNESSES: [(u64, &[u64]); NUM_BASES] = [
        (2_046,                             &[2]                                                        ),
        (1_373_652,                         &[2, 3]                                                     ),
        (9_080_190,                         &[31, 73]                                                   ),
        (25_326_000,                        &[2, 3, 5]                                                  ),
        (4_759_123_140,                     &[2, 7, 61]                                                 ),
        (1_112_004_669_632,                 &[2, 13, 23, 1_662_803]                                     ),
        (2_152_302_898_746,                 &[2, 3, 5, 7, 11]                                           ),
        (3_474_749_660_382,                 &[2, 3, 5, 7, 11, 13]                                       ),
        (341_550_071_728_320,               &[2, 3, 5, 7, 11, 13, 17]                                   ),
        (3_825_123_056_546_413_050,         &[2, 3, 5, 7, 11, 13, 17, 19, 23]                           ),
        (u64::MAX,                          &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]               ),
    ];

    if n == 2 || n == 3 {
        return true;
    } else if n <= 1 || n % 2 == 0 || n % 3 == 0 {
        return false;
    }

    let mut i = 5;
    let l = n.ilog2() as u64;
    while i <= l {
        if n % i == 0 || n % (i + 2) == 0 {
            return false;
        }
        i += 6;
    }

    // 找到 r 使得 n = 2^d * r + 1 对于某些 r >= 1
    let mut d = n - 1;
    while d % 2 == 0 {
        d >>= 1;
    }

    let mut i = 0;
    while i < NUM_BASES && WITNESSES[i].0 < n {
        i += 1;
    }
    let witnesses = WITNESSES[i].1;

    let mut i = 0;
    while i < witnesses.len() && witnesses[i] < n {
        if !miller_test(d, n, witnesses[i]) {
            return false;
        }
        i += 1;
    }

    true
}

#[inline(always)]
const fn miller_test(mut d: u64, n: u64, k: u64) -> bool {
    let mut x = mod_pow(k, d, n);
    if x == 1 || x == n - 1 {
        return true;
    }

    while d != n - 1 {
        x = mod_mul(x, x, n);
        d *= 2;

        if x == 1 {
            return false;
        } else if x == n - 1 {
            return true;
        }
    }

    false
}

/// 计算 (`base` ^ `exp`) mod `modulo` 而不会溢出。
#[inline(always)]
pub const fn mod_pow(mut base: u64, mut exp: u64, modulo: u64) -> u64 {
    let mut res = 1;

    base %= modulo;

    while exp > 0 {
        if exp % 2 == 1 {
            res = mod_mul(res, base, modulo);
        }
        base = mod_mul(base, base, modulo);
        exp >>= 1;
    }

    res
}

/// 计算 (`a` * `b`) mod `modulo` 而不会溢出。
#[inline(always)]
pub const fn mod_mul(a: u64, b: u64, modulo: u64) -> u64 {
    ((a as u128 * b as u128) % modulo as u128) as u64
}
